* Rep_Survival 11 years test of standard errors.sps.
* Produce file of survival times for teeth filled in period.
* Adjusted for time dependency.
* Test standard errors by using random censoring as alternative to weighted adjustment.
* Written by PSKL on 24/09/02.

* Work from file of observed restoration lives.

get file='D:\Longitudinal Data\Rep_Observed Restoration Lives 1991 to 2001.sav'.

* Restrict to restorations placed before the end of the observation period 31/12/2001.

select if (doplacn<yrmoda(2001,12,31)-yrmoda(1899,12,31)).

* Check whether time is in acceptable range, and set reint flag.

compute reint=0.
if not (missing(time)) reint=1.
if (time le 0) reint=0.
if (doplacn+time>yrmoda(2001,12,31)-yrmoda(1899,12,31)) reint=0.

execute.

* Create weights to adjust for time dependent censoring.

* Double up the cases with no re-intervention and weight by probability of re-intervention.
compute caseid=$casenum.
execute.

* Look up probability of reattendance after last visit.
sort cases by interval.
match files file=*
 /table='D:\Longitudinal Data\Rep_reattendance probability.sav'
 /by interval.
compute recid=1.
sort cases by caseid.
save outfile='d:\temp1.sav'.

select if (reint=0).
compute recid=2.
sort cases by caseid.
save outfile='d:\temp2.sav'.

add files file='d:\temp1.sav'
 /file='d:\temp2.sav'
 /by caseid .

Compute weight = 1000.
Do if (reint=0).
if (recid=1) weight = 1000-trunc(1000*prreatt+0.5).
if (recid=2) weight = trunc(1000*prreatt+0.5).
end if.

save outfile='d:\temp3.sav'.


get file='d:\temp3.sav'.

***************************************************.
* Compare and calculate survival function and statistics, and standard errors.


* Assume censored cases are lost at one year (365 days) after last visit.

* recalculate time if reint=0.
* Impute effective time of censoring.
Do if (reint=0).
If (recid=1) time=min(lastvn+365,yrmoda(2001,12,31)-yrmoda(1899,12,31))- doplacn.
if (recid=2) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.
end if.

save outfile= 'd:\temp4.sav'.
get file='d:\temp4.sav'.

*  Test 1 Adjusted KM, but increased weights give incorrect confidence intervals.

weight by weight.

KM
  time  /STATUS=reint(1)
    /PRINT MEAN
    /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
    /Save survival (survadj).

Kaplan-Meier

**Notes**
Output Created		28-MAY-2003 12:18:13
Comments
Input	Data	d:\temp4.sav
	Filter	<none>
	Weight	WEIGHT
	Split File	<none>
	N of Rows in Working Data File	814068
Syntax		KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survadj).
Resources	Elapsed Time	0:01:59.25

 Survival Analysis for TIME




          Survival Time    Standard Error   95% Confidence Interval

 Mean:      2553.18                 .09     (  2553.01,   2553.36 )
 (Limited to   4017.0 )
 Median:    3061.00                 .46     (  3060.09,   3061.91 )



                                                               Percentiles

                     45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

 Value             3808.00      3061.00      2486.00      1997.00      1593.00      1255.00       967.00       723.00       515.00
 Standard Error        .            .46          .33          .25          .19          .15          .12          .09          .07

                                                               Percentiles

                     90.00        95.00

 Value              332.00       175.00
 Standard Error        .05          .

>Warning # 3211
>On at least one case, the value of the weight variable was zero, negative,
>or missing.  Such cases are invisible to statistical procedures and graphs
>which need positively weighted cases, but remain on the file and are
>processed by non-statistical facilities such as LIST and SAVE.


* Test 2 Unadjusted KM, for confidence intervals.

weight off.

* Censor all non-returners.
select if (recid=1).

KM
  time  /STATUS=reint(1)
   /PRINT MEAN
   /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
   /Save survival (survunad).

Kaplan-Meier

**Notes**
Output Created		28-MAY-2003 12:20:18
Comments
Input	Data	d:\temp4.sav
	Filter	<none>
	Weight	<none>
	Split File	<none>
	N of Rows in Working Data File	485400
Syntax		KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survunad).
Resources	Elapsed Time	0:01:32.18




 Survival Analysis for TIME




          Survival Time    Standard Error   95% Confidence Interval

 Mean:      2498.92                2.89     (  2493.26,   2504.58 )
 (Limited to   4017.0 )
 Median:    2824.00               12.00     (  2800.48,   2847.52 )



                                                               Percentiles

                     45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

 Value             3416.00      2824.00      2320.00      1890.00      1519.00      1211.00       941.00       710.00       509.00
 Standard Error        .          12.00         9.20         7.24         5.66         4.44         3.50         2.70         2.04

                                                               Percentiles

                     90.00        95.00

 Value              332.00       175.00
 Standard Error       1.47          .

* Test 3 Censoring using randomisation instead of weighting.

* Censor the proportion prreatt of all non-returners at the end of the observation period (instead of earlier).

Set seed = 10102002.
compute random=uniform(1).
if (reint=0 and random<prreatt) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.

KM
  time  /STATUS=reint(1)
   /PRINT MEAN
   /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
   /Save survival (survrand).

Kaplan-Meier

**Notes**
Output Created		28-MAY-2003 12:21:50
Comments
Input	Data	d:\temp4.sav
	Filter	<none>
	Weight	<none>
	Split File	<none>
	N of Rows in Working Data File	485400
Syntax		KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survrand).
Resources	Elapsed Time	0:03:16.68




 Survival Analysis for TIME




          Survival Time    Standard Error   95% Confidence Interval

 Mean:      2552.64                2.82     (  2547.11,   2558.16 )
 (Limited to   4017.0 )
 Median:    3059.00               14.67     (  3030.25,   3087.75 )



                                                               Percentiles

                     45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

 Value             3807.00      3059.00      2485.00      1996.00      1592.00      1254.00       966.00       723.00       515.00
 Standard Error        .          14.67        10.42         7.99         6.11         4.74         3.67         2.80         2.10

                                                               Percentiles

                     90.00        95.00

 Value              332.00       175.00
 Standard Error       1.49          .

save outfile='D:\temp5.sav'.

* Rerun to test robustness.

get file='d:\temp4.sav'.
select if (recid=1).
compute random=uniform(1).
if (reint=0 and random<prreatt) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.

KM
  time  /STATUS=reint(1)
   /PRINT MEAN
   /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95).

Kaplan-Meier

**Notes**
Output Created		28-MAY-2003 12:26:36
Comments
Input	Data	d:\temp4.sav
	Filter	<none>
	Weight	<none>
	Split File	<none>
	N of Rows in Working Data File	485400
Syntax		KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95).
Resources	Elapsed Time	0:00:43.02




 Survival Analysis for TIME




          Survival Time    Standard Error   95% Confidence Interval

 Mean:      2552.92                2.82     (  2547.40,   2558.44 )
 (Limited to   4017.0 )
 Median:    3059.00               14.67     (  3030.26,   3087.74 )



                                                               Percentiles

                     45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

 Value             3807.00      3059.00      2486.00      1997.00      1593.00      1255.00       967.00       723.00       515.00
 Standard Error        .          14.67        10.42         7.98         6.11         4.75         3.68         2.80         2.10

                                                               Percentiles

                     90.00        95.00

 Value              332.00       175.00
 Standard Error       1.49          .

**************************************************************************************.

* Create an Excel file for comparing survivor functions.
get file='d:\temp5.sav'.

weight off.

select if (reint=1).
sort cases by time.
aggregate outfile=*
 /presorted
 /break time
 / survadj survunad survrand
= first( survadj survunad survrand ).

SAVE TRANSLATE OUTFILE='I:\Research Projects\longevity\Rep_1991 to 2001 adj unadj and rand.xls'
  /TYPE=XLS /MAP /REPLACE /FIELDNAMES.

Data written to I:\Research Projects\longevity\Rep_1991 to 2001 adj unadj and rand.xls.
4 variables and 3775 cases written to range: SPSS.
Variable: TIME       Type: Number   Width:  8   Dec: 2
Variable: SURVADJ    Type: Number   Width: 10   Dec: 5
Variable: SURVUNAD   Type: Number   Width: 10   Dec: 5
Variable: SURVRAND   Type: Number   Width: 10   Dec: 5





* Create an Excel file for calculation of standard errors.

get file='d:\temp4.sav'.

weight by weight.

sort cases by time reint.

save outfile = 'd:\temp6.sav'.
get file='d:\temp6.sav'.

compute censored=0.
if (reint=0) censored=1.

aggregate outfile=*
 /presorted
 /break time
 /reints censored =sum(reint censored ).

>Warning # 3211
>On at least one case, the value of the weight variable was zero, negative,
>or missing.  Such cases are invisible to statistical procedures and graphs
>which need positively weighted cases, but remain on the file and are
>processed by non-statistical facilities such as LIST and SAVE.


compute reints=reints/1000.
compute censored = censored/1000.


SAVE TRANSLATE OUTFILE='I:\Research Projects\longevity\Rep_1991 to 2001 times for se.xls'
  /TYPE=XLS /MAP /REPLACE /FIELDNAMES.

Data written to I:\Research Projects\longevity\Rep_1991 to 2001 times for se.xls.
3 variables and 4010 cases written to range: SPSS.
Variable: TIME       Type: Number   Width:  8   Dec: 2
Variable: REINTS     Type: Number   Width:  8   Dec: 2
Variable: CENSORED   Type: Number   Width:  8   Dec: 2